A WEAK GALERKIN MIXED FINITE ELEMENT METHOD FOR 
BIHARMONIC EQUATIONS 
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Abstract. This article introduces and analyzes a weak Galerkin mixed finite element method 
for solving the biharmonic equation. The weak Galerkin method, first introduced by two of the 
authors (J. Wang and X. Ye) in |52| for second order elliptic problems, is based on the concept of 
discrete weak gradients. The method allows the use of completely discrete finite element functions on 
partitions of arbitrary polygon or polyhedron. In this article, the weak Galerkin method is applied 
to discretize the Ciarlet-Raviart mixed formulation for the biharmonic equation. In particular, an 
a priori error estimation is given for the corresponding finite element approximations. The error 
analysis essentially follows the framework of Babuska, Osborn, and Pitkaranta !8 and uses specially 
designed mesh-dependent norms. The proof is technically tedious due to the discontinuous nature of 
the weak Galerkin finite element functions. Some computational results are presented to demonstrate 
the efficiency of the method. 
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1. Introduction. In this paper, we are concerned with numerical methods for 
the following biharmonic equation with clamped boundary conditions 

A 2 u = f in £1, 
q L j u — on dSl, 

— = on all, 
on 

where fl is a polygonal or polyhedral domain in M. d (d — 2, 3). To solve the problem 
using a primal-based conforming finite element method, one would need C 1 
continuous finite elements, which usually involve large degree of freedoms and hence 
can be computationally expensive. There are alternative numerical methods, for 
example, by using either nonconforming elements [H [351 SI] , the C° discontinuous 
Galerkin method [SH H], or mixed finite element methods [TT1 \W\ |2"01 123 1 133 1 |53j 153 1 
|3H1 [37J 13H SO] ■ One of the earliest mixed formulation proposed for (|1.1[) is the Ciarlet- 
Raviart mixed finite element formulation [50] which decomposes (jl.ll) into a system 
of second order partial differential equations. More precisely, in this formulation, one 
introduces a dual variable w = —Ait and rewrites the four-order biharmonic equation 
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into two coupled second order equations 



(1.2) 



( 



w + Au = 0, 
-Aw = /, 



In [20] , the above system of second order equations is discretized by using the standard 
H 1 conforming elements. However, only sub-optimal error estimates are proved in 
[20j for quadratic or higher order of elements. Improved error estimates have been 
established in [5J 123 EE] f° r quadratic or higher order of elements. In [SJ, Babuska, 
Osborn and Pitkaranta pointed out that a suitable choice of norms are L 2 for w and 
H 2 for u, or equivalent, in order to use the standard LBB stability analysis. In this 
sense, one has "optimal" order of convergence in H 2 norm for u and in L 2 norm for w, 
for quadratic or higher order of elements. However, when equal order approximation 
is used for both u and w, the "optimal" order of error estimate is restricted by the 
interpolation error in H 2 norm, and thus may not be really optimal. Moreover, this 
standard technique does not apply to the piecewise linear discretization, since in this 
case the interpolation error can not even be measured in H 2 norm. A solution to this 
has been proposed by Scholz [55]. Using an L°° argument, Scholz was able to improve 
the convergence rate in L 2 norm for w by hz, and this theoretical result is known to 
be sharp. Also, Scholz's proof works for all equal-order elements including piecewise 
linears. 

The goal of this paper is to propose and analyze a weak Galerkin discretization 
method for the mixed formulation (jl.2[) . The weak Galerkin method was recently 
introduced in [52] for second order elliptic equations. It is an extension of the standard 
Galerkin finite element method where classical derivatives were substituted by weakly 
defined derivatives on functions with discontinuity. Optimal order of a priori error 
estimates has been observed and established for various weak Galerkin discretization 
schemes for second order elliptic equations [SH (53J [32] . A numerical implementation 
of weak Galerkin was discussed in [43] [42] for some model problems. 

Applying the weak Galerkin method to both second-order equations in (|1.2[) ap- 
pears to be trivial and straight-forward at first glance. However, the application 
turns out to be much more complicated than simply combining one weak Galerkin 
scheme with another one. The application is particularly non-trivial in the mathe- 
matical theory on error analysis. In deriving an a priori error estimate, we follow 
the framework as developed in [5J by using mesh-dependent norms. Many commonly 
used properties and inequalities for standard Galerkin finite element method need to 
be re-derived for weak Galerkin methods with respect to the mesh-dependent norms. 
Due to the discrete nature of the weak Galerkin functions, technical difficulties arise 
in the derivation of inequalities or estimates. The technical estimates and tools that 
we have developed in this paper should be essential to the analysis of weak Galerkin 
methods for other type of modeling equations. They should also play an important 
role in future developments of preconditioning techniques for weak Galerkin methods. 
Therefore, we believe this paper provides useful technical tools for future research, in 
addition to introducing an efficient new method for solving biharmonic equations. 

The paper is organized as follows. In Section [2l a weak Galerkin discretization 
scheme for the Ciarlet-Raviart mixed formulation of the biharmonic equation is intro- 
duced and proved to be well-posed. Section [3] is dedicated to defining and analyzing 
several technical tools, including projections, mesh-dependent norms and some esti- 
mates. With the aid of these tools, an error analysis is presented in Section^ Finally, 
in Section [SJ we report some numerical results that show the efficiency of the method. 
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2. A Weak Galerkin Finite Element Scheme. For illustrative purpose, we 
consider only the two-dimensional case of (|1.1[) and the corresponding weak Galerkin 
method will be based on a shape- regular triangulation of the domain ft. 

Let D C fl be a polygon, we use the standard definition of Sobolev spaces H S (D) 
and Hq(D) with s > (e.g., see [1] [21] for details). The associated inner product, 
norm, and semi-norms in H S (D) are denoted by (-,-) s ,d, || • || s ,d, and | • | r ,_D,0 < 
r < s, respectively. When s = 0, H°(D) coincides with the space of square integrable 
functions L 2 (D). In this case, the subscript s is suppressed from the notation of norm, 
semi-norm, and inner products. Furthermore, the subscript D is also suppressed when 
D = fl. For s < 0, the space H S (D) is defined to be the dual of H^ S (D). 

Occasionally, we need to use the more general Sobolev space W s,p (fl), for 1 < 
p < oo, and its norm || ■ ||w>>p(n)- The definition simply follows the standard one given 
in [TJ [H]. When s = 0, the space W s ' p (fl) coincides with LP(fl). 

The above definition/notation can easily be extended to vector- valued and matrix- 
valued functions. The norm, semi-norms, and inner-product for such functions shall 
follow the same naming convention. In addition, all these definitions can be trans- 
ferred from a polygonal domain D to an edge e, a domain with lower dimension. 
Similar notation system will be employed. For example, || • \\ s . e and || • || e would 
denote the norm in H 8 {e) and L 2 (e) etc. We also define the H(div) space as follows 

H(div, fl) = {q : q £ [L 2 (fl)} 2 , V • q £ L 2 (fl)}. 

Using notations defined above, the variational form of the Ciarlet-Raviart mixed 
formulation (|1.2[> seeks u 6 H^fl) and w £ 1 (SI) satisfying 

f( W ,0)-(Vu,V0)=O for all e H^Sl), 

\ (Vw, V-0) = (/, V) for all i> e ffl(n). 

For any solution w and u of (|2.1j) . it is not hard to see that w = —Au. In addition, 
by choosing <f> = 1 in the first equation of (|2.1[) . we obtain 

wdx = 0. 

Define ti^fl) C by 

i? 1 (rj) = {w: Bgtff!)), / vdz^O}, 

Jo 

which is a subspace of _ff 1 (f2) with mean- value free functions. Clearly, the solution w 
of (j2~T|) is a function in ^^fi). 

One important issue in the analysis is the regularity of the solution u and w. 
For two-dimensional polygonal domains, this has been thoroughly discussed in |12j . 
According to their results, the biharmonic equation with clamped boundary condition 
(|l.ip satisfies 

(2-2) HU- fc < e||/||_ fc! 

where c is a constant depending only on the domain Q. Here the parameter k is 
determined by 

k = 1 if all internal angles of f) are less than 180° 

k = if all internal angles of fl are less than 126.283696 • ■ ° 
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The above regularity result indicates that the solution u G H 3 (Q) when fl is a convex 
polygon and / G — 1 (f2). It follows that the auxiliary variable w G Moreover, 
if all internal angles of Q are less than 126.283696 ■ • ° and / G L 2 (n), then u G ff 4 (ft) 
and to G H 2 (Vt). The drawback of the mixed formulation (|2.1I) is that the auxiliary 
variable u> may not possess the required regularity when the domain is non-convex. 
We shall explore other weak Galerkin methods to deal with such cases. 

Next, we present the weak Galerkin discretization of the Ciarlet-Raviart mixed 
formulation. Let 7ft be a shape-regular, quasi-uniform triangular mesh on a polygonal 
domain fi, with characteristic mesh size h. For each triangle K G 7ft, denote by Kq 
and dK the interior and the boundary of K, respectively. Also denote by Kk the 
size of the element K. The boundary dK consists of thee edges. Denote by Eh the 
collection of all edges in 7ft. For simplicity of notation, throughout the paper, we use 
"<" to denote "less than or equal to up to a general constant independent of the mesh 
size or functions appearing in the inequality" . 

Let j be a non-negative integer. On each K G 7ft, denote by Pj(Ko) the set of 
polynomials with degree less than or equal to j. Likewise, on each e£4, Pj{ e ) is the 
set of polynomials of degree no more than j . Following [52] , we define a weak discrete 
space on mesh Th by 

Vft = {v : v\ Ko G Pj(K ), K G v\ e G Pj (e), e G S h }. 

Observe that the definition of Vh does not require any continuity of w £ 14 across 
the interior edges. A function in Vh is characterized by its value on the interior of 
each element plus its value on the edges/faces. Therefore, it is convenient to represent 
functions in Vh with two components, v = {vo,Vb}, where vq denotes the value of v 
on all Kq and Vj, denotes the value of v on £h ■ 

We further define an L 2 projection from onto Vft by setting QhV = 

{QqV, Qbv}, where Qov\k is the local L 2 projection of v in Pj(K ), for K G 7ft, 
and Qi,v\ e is the local L 2 projection in Pj(e), for e G £h- To take care of the homoge- 
neous Dirichlet boundary condition, define 

Vb,ft = {v G Vft : v = on £ h n dfl}. 

It is not hard to see that the L 2 projection Qh maps Hq(VL) onto Vb,ft. 

The weak Galerkin method seeks an approximate solution [uh] Wh] G Vo,ft X Vh 
to the mixed form of the biharmonic problem (jl.2l) . To this end, we first introduce a 
discrete L 2 -equivalent inner-product and a discrete gradient operator on Vft. For any 
Vh = {^Oj^Ei} and 4>h = {<j>o,(pb} in Vft, define an inner-product as follows 

{vh,<ph))= ^ ( w o,0o)k+ h K (v Q - Vb,(po - <i>b)dK- 
KeT h Ker h 

It is not hard to see that ((vh, Uft)) = implies Vh = 0. Hence, the inner-product is well- 
defined. Notice that the inner-product ((-, •)) is also well-defined for any v G 
for which v — v and Vb\ e = v\ e is the trace of v on the edge e. In this case, the 
inner-product ((•,•)) i s identical to the standard L 2 inner-product. 

The discrete gradient operator is defined element- wise on each K G 7ft. To this 
end, let RTj(K) be a space of Raviart-Thomas element [U] of order j on triangle K. 
That is, 



RT j (K) = (P j (K)) 2 +xP j (K). 
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The degrees of freedom of RTj (K) consist of moments of normal components on each 
edge of K up to order j, plus all the moments in the triangle K up to order (j — 1). 
Define 

= {qe (L 2 (n)) 2 : q\ K E RT,(K), K E %}■ 

Note that E/j is not necessarily a subspace of H(div, ft), since it does not require any 
continuity in the normal direction across any edge. A discrete weak gradient [52] of 
v h = {vQ^Vb} £ Vh is defined to be a function V w Vh E Eft such that on each K E 7fe, 

(2.3) (V d v h ,q) K = -(vo,V-q)K + (vb,q.-n) dK , for all q E RTj (K) , 

where n is the unit outward normal on dK. Clearly, such a discrete weak gradient is 
always well-defined. Also, the discrete weak gradient is a good approximation to the 
classical gradient, as demonstrated in [52] : 

Lemma 2.1. For any Vh = {vq, Vb} E Vh and K E 7ft, V«,Uh|ir — if and only if 
vq = Vb = constant on K . Furthermore, for any v E H m+1 (Q) , where < m < j + 1, 
we /icwe 

HV^CQhv)- V«|| </i ro |H| m+x . 



We are now in a position to present the weak Galcrkin finite element formulation 
for the biharmonic problem (jl.2p in the mixed form: Find Uh — {ito, u&} E Vq^ and 
w h — {u>o, uib} E Vh such that 

(2 4) )i Wh -> M - (V w Ufc, V w ^>h) = 0, for all ^ = {0 O , 4>b] E 14, 

^(VtoiUh, V w t/Jh) = {f, tpo), for all ip h = {ip , V'fc} € Vox 



Theorem 2.2. The weak Galerkin finite element formulation has one and 

only one solution [uh',Wh] in the corresponding finite element spaces. 

Proof. For the discrete problem arising from (|2.4[) , it suffices to show that the 
solution to (|2.4p is trivial if / = 0; the existence of solution stems from its uniqueness. 

Assume that / = in (|2.4[) . By taking <\>h — Wh and iph = Uh m (|2.4I) and 
adding the two resulting equations together, we immediately have ([wh, WhJ) = 0, 
which implies Wh = 0. Next, by setting <f>h — Uh in the first equation of (|2.4|) . we 
arrive at (V w Uh, VwUh) = 0. By using Lemma |2~T| we see that Uh must be a constant 
in fi, which together with the fact that Uh = on dil implies Uh = in O. This 
completes the proof of the theorem. □ 

One important observation of (|2.4[) is that the solution Wh has mean value zero 
over the domain Q, which is a property that the exact solution w = —Au must possess. 
This can be seen by setting <fih = 1 in the first equation of (|2.4p . yielding 

(w h , 1) = (w h , 1)) = (V w u h , V w l) = 0, 

where we have used the definition of ((•,•)) an d Lemma 12.11 For convenience, we 
introduce a space Vh C Vh defined as follows 

V h = {v h : v h = {v , v b } E Vh, / v dx = 0}. 
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3. Technical Tools: Projections, Mesh-dependent Norms and Some Es- 
timates. The goal of this section is to establish some technical results useful for 
deriving an error estimate for the weak Galerkin finite element method (|2.4I) . 

3.1. Some Projection Operators and Their Properties. Let P/j be the 

L 2 projection from (L 2 (f2)) 2 to S/ l7 and 11^ be the classical interpolation [16] from 
(iJ 7 (17)) 2 ,7 > |, to Efj defined by using the degrees of freedom of Ti h in the usual 
mixed finite element method. It follows from the definition of 11^ that Il^q G 
H(div,&V) n Eh for all q G (_ff 7 (fi)) 2 . In other words, n^q has continuous nor- 
mal components across internal edges. It is also well-known that II;, preserves the 
boundary condition q ■ n|gji = 0, if it were imposed on q. The properties of II/j 
has been well-developed in the context of mixed finite element methods [TBI EO] . For 
example, for all q G (W m ' p (fl)) 2 where | < m < j + 1 and 2 < p < oo, we have 

(3.1) Qo(V • q) = V • n h q, if in addition q G H(div, fi), 

(3.2) ||q - U h q\\ L p(n) < h m \\q\\ wm , P{n) . 

It is also well-known that for all < m < j + 1 , 
(3-3) llq-PhQll </i m ||q|U- 

Using the above estimates and the triangle inequality, one can easily derive the fol- 
lowing estimate 

(3.4) \\U h Vv-P h Vv\\<h m \\v\\ m+1 
for all v G H m+1 (fl) where \ < m < j + 1. 

Next, we shall present some useful relations for the discrete weak gradient V,„, 
the projection operator P/,, and the interpolation U^. The results can be summarized 
as follows. 

Lemma 3.1. Let 7 > ^ be any real number. The following results hold true, 
(i) For any v G H {ft), we have 

(3.5) V t0 (Q A t;)=Pfc(Vt;). 

(ii) For any q G {H 1 (fl)) 2 fl H(div, Cl) and Vh — {«o, vi,} G Vh, we have 

(3.6) (V • q, W ) = -(n h q, V w v h ) + ^ ((U h q) ■ n,v b } e . 

In particular, if either Vh G Vo t h or q • n = on dQ, then 

(3.7) (V • q, v ) = -(n h q, V w v h ). 

Proof. To prove (|3.5[) . we first recall the following well-known relation [16] 

V • RT^K) = Pj{K ), RT 3 (K) • n| e = P 3 (e). 

Thus, for any w G and K G 7ft, by the definition of \7 W and properties of the L 2 
projection, we have 

(VwQhV, w)k = ~(Qov, V • w) A - + (QbV, w • n) dK 
= -(v, V • w)k + (v, w • n)aK 
= (Vu,w)jc 
= (PfcVu.w)*, 
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which implies (|3.5|) . As to (|3.6|) . using the fact that V-RTj(K) = Pj(K ), the property 
(|3.ip , and the definition of V,,, we obtain 

(V • q, v ) = (Q (V • q), wo) = (V • n, iq , wo) 

= - E ( n '^' ^wVh)K + E ( u 6' n 'iq- n W 

= - E ( n ' lC l' ^wVhlK + E (( n /iq) ' n , Wb) e . 
ifeTh eeT h ndn 

This completes the proof of (|3.6I) . The equality (|3.7I) is a direct consequence of p. 61) 
since the boundary integrals vanish under the given condition. □ 

3.2. Discrete Norms and Inequalities. Let Vh — {vo,Vb} £ Vh- Define on 
each K £ Th 

hh\\l, h , K = \\ v o\\l,K + h \\ v ~ V b \\dK, 
\\Vh\\l,h,K = \\ v o\\l,K + h^Wvo - V b f dK , 
W\\,h, K = Ml,K + ^IN - V b \\ 2 dK . 

Using the above quantities, we define the following discrete norms and semi-norms 
for the finite element space Vh 




W\l.h ■= E \ V h\i,h,K 
\KET h 

It is clear that ||vfi||o/i = Hence, || ■ \\o,h provides a discrete L? norm for 

Vh. It is not hard to see that | • |i h and || • define a discrete H 1 semi-norm and 
a norm for Vh, respectively. Observe that = if and only if Vh = constant. 

Thus, | • is a norm in V) t h and Vh- 

For any K £ Th and e being an edge of K, the following trace inequality is 
well-known 

(3-8) \\9\\l<h-%\& + h*'- 1 \g\l K , \<s<l, 

for all g £ H l (K). Here |^| s ,j<- is the semi-norm in the Sobolev space H S (K). The 
inequality p.8[) can be verified through a scaling argument for the standard Sobolev 
trace inequality in H s with s 6 If g is a polynomial in K, then we have from 

(|3.8[) and the standard inverse inequality that 

(3-9) \\9\\l<h- l \\gf K . 

From (|3.9|) and the triangle inequality, it is not hard to see that for any Vh £ Vh 
one has 
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In the rest of this paper, we shall use the above equivalence without particular men- 
tioning or referencing. 

The following Lemma establishes an equivalence between the two semi-norms 
| • \ 1>h and \\V W ■ ||. 

Lemma 3.2. For any Vh = {v$,v b } G Vh> we have 
(3-10) \v h \i,h < HV^fcH < \v h \i,h- 



Proof. Using the dehnition of \7 W , integration by parts, the Schwarz inequality, 
the inequality (|3.9[) . and the Young's inequality, we have 

llVtuW/tHl- = -(vq, V • V w v h ) K + (v b ,V w v h ■ n) dK 
= {v b - v , V w v h ■ n) dK + (Vv Q ,V w v h ) K 

< ho - Vb\\dK\\^wVh ■ n\\sK + ||V«o||if ||V w Ufc||K 

< \\v - v b \\ dK h~i\\V w v h \\ K + \\Vvo\\ K \\V w v h \\K 

< \\VwVh\\K (||Vw |k + h~i\\v - v b \\ dK ) ■ 

This completes the proof of ||V TO U/i|| < 

To prove < ||V w «/i||, let K G Th be any element and consider the following 

subspace of RTj(K) 

D(J, K) := {q G RTj(K) : q • n = on dK}. 

Note that D(j, K) forms a dual of {K)f . Thus, for any Vv G (Pj-^K)) 2 , one 
has 

(3.11) ||V«o||k = sup n I, ■ 

q£.D(j,JC) iRllK 

It follows from the integration by parts and the dehnition of \7 W that 

(Vv ,q) K = -(v , V- q)if = (V^Vfc, q)*r, 
which, together with p. Ill) and the Cauchy- Schwarz inequality, gives 

(3.12) ||Vw |k < ||V„w/i||jr. 

Note that for j = 0, we have Vvq = and the above inequality is satisfied trivially. 

Analogously, let e be an edge of K and denote by D e (j,K) the collection of all 
q G RTj(K) such that all degrees of freedom, except those for q • n| e , vanish. It is 
well-known that D e (j, K) forms a dual of Pj(e). Thus, we have 

, q1 q\ ii ii (w -Wfc,q-n) e 

(3.13) ||«o-«b||e= sup 



qeD c (],K) l|q ,n l|e 

It follows from (|2.3|) and the integration by parts on (vq, V • q}x that 

(3.14) (V w v h ,q) K = (Vv Q ,q) K + (v b - V ,q.- n) aK , V q G RTj(K). 
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In particular, for q £ D e (j,K), we have 

(V«o, q) K = 0, (v b - v 0l q ■ n)o K = (vb - i>o,q • n) e . 
Substituting the above into Q3.14p yields 

(3.15) (VwVh,q)K = {vb -«o,q- n) e , V q G D e (j,K). 

Using the Cauchy-Schwarz inequality we arrive at 

\{vb -u ,q-n) e | < \\\7 w v h \\ K \\q\\ K , 

for all q £ D e (j,K). By the scaling argument, for such q £ D e (j,K), we have 
WiWk S 1| q ■ n|| e . Thus, we obtain 

\(v b - v Ql q- n) e \ < h?\\V w v h \\ K ||q- n|| e , Vq £ D e (j,K), 

which, together with (|3.13l) . implies the following estimate 

\\ v - v b\\e < h^\\V w Vh\\K- 

Combining the above estimate with f|3 . 1 2[) gives a proof of ^ < HV^w^H. This 
completes the proof of (I3.10|) . □ 

The discrete semi-norms satisfy the usual inverse inequality, as stated in the 
following Lemma. 

Lemma 3.3. For any Vh = {vo, Vb} £ Vh, we have 
(3-16) \v h \ ljh < /»- 1 ||«fc||o,fc. 

Consequently, by combining i3. 1 0\) and i3.16]) . we have 

(3.17) ||V„Wfc||<fe- 1 |K||o,fc. 

Proof. The proof follows from the standard inverse inequality and the definition 
of || • | [ o,/* and | • |i^; details are thus omitted. □ 

Next, let us show that the discrete semi- norm ||V lu (-)||, which is equivalent to 
• \\.h as proved in Lemma 13.21 satisfies a Poincare-type inequality. 

Lemma 3.4. The Poincare-type inequality holds true for functions in Vb,/i and 
Vh- In other words, we have the following estimates: 

(3.18) \\v h \\ , h <\\V w v h \\ Vv h £V , h , 

(3.19) \\v h \\ , h < \\V w v h \\ Vv h £V h . 

Proof. For any Vh £ Vo,h, let q £ (iJ 1 (f2)) 2 be such that V • q = v and ||q||i < 
|| «o II- Such a vector- valued function q exists on any polygonal domain [3]. One 
way to prove the existence of q is as follows. First, one extends Vh by zero to a 
convex domain which contains Q. Secondly, one considers the Poisson equation on 
the enlarged domain and set q to be the flux. The required properties of q follow 
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immediately from the full regularity of the Poisson equation on convex domains. By 
(|3.1[) . we have 

llHfcqH < ||q||i < ||«o||. 
Consequently, by (13.71) and the Schwarz inequality, 

boll 2 = («b,V-q) = -(n h q,V„» h ) < ||«b||||V„w h ||. 
It follows from Lemma T3. 21 that 

Y h W v ° ~ v b\\dK < Y h ~ l \\vo-i>b\\aK < \vh\l, h < ||V w « h || 2 . 
Ken Ken 

Combining the above two estimates gives a proof of the inequality (|3.18|> . 

As to (|3.19p . since Vh 6 Vh has mean value zero, one may find a vector- valued 
function q satisfying V • q = uo and q • n = on dft (see [3] for details) . In addition, 
we have ||q||i < \\vo\\- The rest of the proof follows the same avenue as the proof of 
HUSK. □ 

Next, we shall introduce a discrete norm in the finite element space Vo t h that 
plays the role of the standard H 2 norm. To this end, for any internal edge e £ £h, 
denote by K\ and K% the two triangles sharing e, and by ni, 112 the outward normals 
with respect to K\ and K%- Define the jump on e by 

{V w iph ■ n] = {V w ip h )\ Kl ■ ni + (V w iph)\K 2 ■ n 2 . 

If the edge e is on the boundary <9f2, then there is only one triangle K which admits 
e as an edge. The jump is then modified as 

[V«i^ • n] = (y w iph)\K ■ n. 

For tph€Vo jh , define 

(3.20) ll^lll = ( Y IIV-V W ^||^+ Y ^IllV^-n]^) 

\K£T h ee£ h 




Lemma 3.5. The map \\\ ■ \\\ : Vo.h — > R, as given in I13.20\) , defines a norm in the 
finite element space Vo,h- Moreover, one has 

(3.21) (V^fc, V w i/> h ) < \\v h \\^ h \\^ h \\ V v h 6 V h , fa £ V ,h, 

(3.22) sup n — n > \\\^h\\\ V ^ft e Vo,h. 

v h ev h Ph\\o t h 



Proof. To verify that ||| • ||| defines a norm, it is sufficient to show that \\\iph\\\ = 
implies iph = 0. To this end, let \\\iph\\\ = 0. It follows that V • V w iph = on each 
element and [V„,V/t" n ] = on each edge. The definition of the discrete weak gradient 
Vuj then implies the following 

(V w ip h , V^Vfo) = Y (~(^o, V • V w 4)/f + (^ fc , V w ip h ■ n)aAr) = 0. 
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Thus, we have V^^/i = 0. Since ip^ £ Vo,h, then V w iph — implies iph = 0. This 
shows that ||| • ||| defines a norm in Vo.h- The inequality (|3.21[) follows immediately 
from the following identity 

(V w v h , V w tph) = Y (-(^o, V • V w iph)K + (vb, 'Vwiph ■ n)ax) 
Ken 

and the Schwarz inequality. 

To verify (|3.22p . we chose a particular u? S such that 

«o = -V • V u) i/'/i in K , 

v b = h^lVwiph ■ n] on edge e. 

It is not hard to see that ||?^llo,fc, ^5 IIIV^II- Thus, we have 

v h ev h \\vh\\o,h \K\\o,h 

_ Y<Ken (~( u 0' v ' v ^h)if + (uj, V^Vfc ' n W) 



Kllo.h 



ll^ll 2 I,, , I, 

£ III II 



This completes the proof of the lemma. □ 

Remark 3.1. Using the boundedness !i3.21\) and the discrete Poincare inequality 
k3. 1 8\) we have the following estimate for all iph G Vo,h 

llv^ll 2 = {Vwi>h,Vwi>h) < Hh\\o.h\\\iPh\\\ < || v w i/> h II Uh |||. 

This implies that \[V w iph\\ < \\\iph\\\ ■ In other words, \\\ ■ \\\ is a norm that is stronger 
than || • ||i, ft. In fact, the norm \\\ ■ \\\ can be viewed as a discrete equivalence of the 
standard H 2 norm for smooth functions with proper boundary conditions. 

Next, we shall establish an estimate for the L 2 projection operator Qh in the 
discrete norm || • \\o,h- 

Lemma 3.6. Let Qh be the L 2 projection operator into the finite element space 
Vh- Then, for any v G H m (£l) with i < m < j + 1, we have 

(3.23) \\v-Q h vh,h<h m \\v\\ m . 



Proof. For the I? projection on each element K, it is known that the following 
estimate holds true 

(3.24) \\v-Q v\\ K <h m \\v\\ m . K . 

Thus, it suffices to deal with the terms associated with the edges/faces given by 

(3.25) Yl h W( v ~ - ( u ~ Qbv)\\ 2 dK = h WQ° v ~ Qbv\\ 2 dK . 

K K 

Since Qb is the L 2 projection on edges, then we have 

\\Q V - Qbv\\l K < \\ v ~ Qov\\dK- 
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Let s G (i, 1] be any real number satisfying s < m. It follows from the above inequality 
and the trace inequality (|3.8[) that 

WQav - Qbv\\ 2 9K < h^Wv - Q v\\ 2 K + h 2 *- 1 ^ - Qo< K - 
Substituting the above into (I3.25|) yields 

Y,h\\{v - Q Q v) -(v- Q b v)\\ 2 dK < £ - Qov\\% + h 2s \v - Q v\l K ) 

K K 

<^ 2m IMlL 

which, together with (|3.24p . completes the proof of the lemma. □ 

3.3. Ritz and Neumann Projections. To establish an error analysis in the 
forthcoming section, we shall introduce and analyze two additional projection opera- 
tors, the Ritz projection Rh and the Neumann projection Nh, by applying the weak 
Galerkin method to the Poisson equation with various boundary conditions. 

For any v G H^(fi) n i/ 1+7 (Sl) with 7 > i, the Ritz projection R h v G V . h is 
defined as the unique solution of the following problem: 

(3.26) (V w (R h v),V w ^ h ) = (U h Vv,V w ip h ), V e V , h . 

Here 7 > 4 in the definition of Rh is imposed to ensure that II/jVu is well-defined. 
From the identity (|3.7[) . clearly if Av € L 2 (f2), then RhV is identical to the weak 
Galerkin finite element solution 52 to the Poisson equation with homogeneous Dirich- 
let boundary condition for which v is the exact solution. Analogously, for any 
v G H 1 ^) n i/ 1+7 (f2) with 7 > 5, we define the Neumann projection N^v G Vh 
as the solution to the following problem 

(3.27) (V w {N h v), V w Tp h ) = (IlftVu, Vt,,^), V^e %. 

It is useful to note that the above equation holds true for all iph G Vh as V^l = 0. 
Similarly, if Av G L 2 (f2) and in addition dv/dn — on dfl, then NhV is identical to 
the weak Galerkin finite element solution to the Poisson equation with homogeneous 
Neumann boundary condition, for which v is the exact solution. The well-posedness of 
Rh and Nh follows immediately from the Poincare-type inequalities (I3.18|) and (|3.19[) . 
Using (|3.5p . it is easy to see that for all iph S Vo,h we have 

(3.28) (V w (Q h v - Rhv), V w ip h ) = ((Pj, - n fc )V«, V w ip h ). 
And similarly, for all iph £ Vh, 

(3.29) (V w (Qh,v - N h v), V w Tp h ) = ((Ph - U h )Vv, V w ip h ). 

From the definitions of Vh and Qh, clearly Qh maps _ff 1 (f2) into Vh- 
For convenience, let us adopt the following notation 

{R v, R b v} := R h v, {N Q v, N b v} := N h v, 

where again the subscript "0" denotes the function value in the interior of triangles, 
while "6" denotes the trace on £h- For Ritz and Neumann projections, the following 
approximation error estimates hold true. 
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Lemma 3.7. For v e H%(Sl) n H m+1 (Q) or H^fl) n ff m+1 (0), w/iere ± < m < 
j + 1. we have 

(3.30) ||V w (Q h t; - R h v)\\ < h m \\v\\ m+1 , 

(3.31) ||V„(Q fc t;- JV fc «)[| <h m \\v\\ m+1 . 

Moreover, assume Av G L 2 (17) cme? i/iai </ie Poisson problem in £1 uwf/i either the 
homogeneous Dirichlet boundary condition or the homogeneous Neumann boundary 
condition has H 1+s regularity, where h < s < 1, i/ien 

(3.32) ||Q„« - i? «|| < h m+s \\v\\ m+1 + h 1+s \\(I - Q )A«||, 

(3.33) ||Qo« - N v\\ < h m+min ^+^\\v\\ m+1 + h 1+s \\(I - Q )Av\\. 

Proof. The estimates (|330>(|33Tj) follow immediately from (|3T28]) - ((3T29|) . ([3~4)l . 
and the Schwarz inequality. Next, we prove p.33|) by using the standard duality 
argument. Let <fi G H 1 ^) be the solution of — A<fi — Qqv — Nqv with boundary 
condition J-j =0. Note that 4> is well-defined since QhV — NhV £ Vh- According to 

the regularity assumption, we have 4> £ H 1+S (£l) and ||</>||i+ s < \\Qqv — Nqv\\. Then, 
by (|3.7p . ()3.29|) . the Schwarz inequality and (|3.4j) . we arrive at 

||Qo« - N v f = (Q v - N Q v, -A<£) = (n fc V<£, V w (Q fc v - N h v)) 

= (II h V<t> - V w {N h <l>), V w {Q h v - N h v)) + ((P h - U h )Wv, V w (N h <f>)) 



< (J|n fc V0 - P h V^|| + \\V w (Qh<f> - AW) || J ||V W (Q^ - A^)|| 

+ ((P h - IL h )Vv, V w {N h <t> - Q h <t>)) + {(P h - IL h )Vv, P h V0) 

< h m+s \\<j>\\i +a \Mm+i + ((/ - n h )Vv, P h V<f>). 

Using integration by parts, the triangular inequality and the definition of 11^, we have 

((j-n^Vu.PftV^) 

=((J - U h )Vv, (P h - I)V<f>) + ((I - U h )Vv, v<t>) 
<h m+s U\\ 1+s \\vln+i + {{I - U h )Vv • n, 0) 9n - (V • (J - Tl h )Vv, 0) 
(3.34) =^ m+s ||0||i+ s ||«|| m +i + ((I - Tl h )Vv .n,4>- Q b <f>) 9 n - ((I - Qo)Av, 0) 
<^+ s ||^||i+ s ||^|U+i + (^-*II^IU+i,ao)(^ millCs+ '' J ' +1) II^IL+i,ao) 

-((I-Qo)Av,(I-Q )<i>) 
</ l ^-daW+J)||0|| 1+ ,|| t ,|| Tn+1 + / l i+-||^|| 1+a ||(j_Q o )A«||. 

In the proof of p.34j) . we have used the fact that Il^Wn) is exactly the L? projection 
of Vt> • n on dQ. Combining the above gives 

||Qo» - N v\\ 2 < (h m+mi < s ^\\v\\ n+1 + h 1+s \\(I - Q )Av\^j U\\ 1+s 

< ^ m+m in( SJ+ i)|| w || m+i + h l+.^j _ Q ) Av \\^ \\Q oV _ NqV \1 

This completes the proof of the estimate (|3.33p . The inequality (I3.32p can be verified in 
a similar way by considering a function <f> 6 Hq(Q,) satisfying a Poisson equation with 
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homogeneous Dirichlet boundary condition. Observe that in this case, the boundary 
integral ((/ — II^)Vu • n, <p)ga in inequality (I3.34j) shall vanish due to the vanishing 
value of 4>. □ 

Remark 3.2. It is not hard to see from \3.3^ that for the Neumann projection, 
if in addition we have §^ — on <9ft, then the term ((I — II/j)Vi; ■ n, <p)gn vanishes 
and one obtains the optimal order estimate of h m+s instead of /j m + min ( s >-J+3) for the 
Neumann projection operator. 

Remark 3.3. If the Poisson equation has the full if 2 regularity in ft, then for v 
satisfying the assumptions of Lemma \3. 7[ we have 



WQqv - R v\\ < h m+1 \\v\\ m+1 + h 2 \\(I - Q )Av\\ for±<m<j + l, 

,| _ N I, < (h m +i\\v\\ m+1 + h 2 \\(I-Q )Av\\ forj = 0, § <m< 1, 

° ll ~\h m + 1 \\v\\ m+1 + h 2 \\(I-Q )Av\\ forj>l,i<m<j+l. 

Again, if in addition, — on dQ, then the Neumann projection has optimal order 
of error estimates, even for j = 0. 

Remark 3.4. The duality argument used in Lemma \3.7\ works only for \\Qqv — 
Rqv\\ and \\Qqv — Nqv\\ . For \\Qh,v — Rh,v\\a,h an d \\Qh,v — Nhv\\o,h involving element 
boundary information, we currently have only sub-optimal estimates. More precisely, 
for v satisfying the assumptions in Lemma \3. 7\ the following estimates hold true. 

\\Q h v - R h v\\ , h < \\V w (Q h v - R h v)\\ < h m \\v\\ m+1 for-<m<j + l, 
(3.35) \ 

WQhv - N h v\\ , h < \\V w (Q h v - N h v)\\ < h m \\v\\ m+1 for - < m < j + 1. 

Although numerical experiments in suggest an optimal order of convergence in 
the || • ||o./i norm, it remains to see if optimal order error estimates hold true or not 
theoretically. 

Another important observation is that, for sufficiently smooth v, V ' w RhV is iden- 
tical to the mixed finite element approximation of Vi>, discretized by using RTj and 
discrete Pj elements. Indeed, we have the following lemma: 

Lemma 3.8. For any v £ Hi n if 1+7 (ft) with 7 > \ and Av £ L 2 (fl), let 
q/i £ Eh n H(div, ft) and Vq £ L 2 (ft) be piecewise Pj polynomials solving 



(3.36) 



{m,x h ) - (v • x/i, «o) = Vxh e E*. n H(div, ft), 

(V ■ q^, ipo) = (Av,ipo) V^o S L 2 (fl) piecewise Pj polynomials. 



In other words, q/j and vq are the mixed finite element solution, discretized using the 
RTj element, to the Poisson equation with homogeneous Dirichlet boundary condition 
for which v is the exact solution. Then, one has V w RhV = q^. 

Proof. We first show that V w Ruv £ E^ n H(div, ft) by verifying that (V w Rhv) ■ n 
is continuous across internal edges. Let e £ £h\(9ft be an internal edge and K\, K2 
be two triangles sharing e. Denote ni and 112 the outward normal vectors on e, with 
respect to K\ and K2, respectively. Let iph € Vo,h satisfy ipt,\e 7^ and ipg, ipi, vanish 
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elsewhere. By the definition of Rh, V ffi and the fact that n^Vi) G H(div, Q), we have 

= (U h Wv - V w R h v, V„Vh) 
= (U h Vv - V w R h v, W w iph)Ki + (n^Vv - V w R h v, W w iph)K 2 
= ((IlfcVtt - VwRhV^Ki ■ ni + (II^Vu - V w -Rftu)|ftr 2 • n 2 ,^ b ) e 
= -(V w Rhv\K 1 ■ ni + V w Rhv\K 2 • n 2 ,ipb)e- 

The above equation holds true for all ipb\ e G Pj(e). Since V^-R^vl^ •ni + V to i?ft,u|K 2 • 
n2 is also in Pj(e), therefore it must be 0. This completes the proof of \7 w RhV G 
H(div,n). 

Next, we prove that V w RhV is identical to the solution q/j of p.36[) . Since the 
solution to (|3.36[) is unique, we only need to show that V^-R/jW, together with a certain 
vq, satisfies both equations in p.36[) . Consider the test function G Vo,h with the 
form tph — {-0o,O}. By the definition of V w , equations (|3.26p and (I3.7p . we have 

(V • V w R h v, ^o) = -{V w R h v, V w <4> h ) = -(U h Vv, V m i) = (Av, Vo)- 

Hence V w Rtv satisfies the second equation of (|3.36p . Now, note that V- is an onto 
operator from n H(div, f2) to the space of piecewise Pj polynomials, which allows 
us to define a vq that satisfies the first equation in (|3.36[) with set to be V w Rh,v. 
This completes the proof the the lemma. □ 

Remark 3.5. Using the same argument and noticing that ^3.27ty holds for all 
?Ph € Vh, one can analogously prove that for v G H 1 ^) D H 1+ ~<(n) with 7 > \ and 

Av G L 2 (n), 

y w N h v G E h H H(div, Q), 

and 

V • V w N h v = Q Av. 

Because V ' w RhV is identical to the mixed finite element solution to the Poisson 
equation, by [50, 30] . we have the following quasi-optimal order L°° estimate: 

(3.37) \\Vv -V w R h v\\ La , (n) <h n+1 \\nh\\\Av\\ Wn ^ {n)l 

for < n < j. Furthermore, for j > 1 and v G W^ +2: °° (fl), we have the following 
optimal order error estimate 

(3.38) ||V« - V w Rhv\\ L ~ m < h n+1 \\v\\ wn+2 ^ m , 
for 1 < n < j. 

Inspired by [15], using the above L°° estimates we obtain the following lemma, 
which will play an essential role in the error analysis to be given in the next section. 

Lemma 3.9. The following quasi-optimal and optimal order error estimates hold 
true: 

(1) Let0<n<j and v G # 1 (ft)rW l+2 < oo (f7). Then for all (j) h = {v ,v b } G V h , 
we have 

(3.39) |(n fc Vv- Vv,RhV,V w <M\ < /i n+ 5|ln/ l ||| ? ;|| w „ + 2,oo (o) ||0 /l |t o . /j . 
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(it) Let j > I, 1 < n < j, and v G ff^O) n W n+2 -°°{tt). Then, for all <f> h 
{vq,vi,} G Vh we have 



(3.40) 



CUkVw 



,RhV,V w (j)h)\ < h n+ 2 \\v\\ Wn +2,co(Q)\\(j>h\\0,h- 



Proof. We first prove part (i). Denote by £,gn the set of all edges in Sh n <9f2. For 
any e G let K e be the only triangle in Th that has e as an edge. Denote by 7an 
the set of all K e , for e G Son- For simplicity of notation, denote q/j = II/iVv — \7 w Rh,v. 
Since (EEftVu — V w RhV,V w tph) = for all V/i G Vo,h, without loss of generality, we 
only need to consider <f>h that vanishes on the interior of all triangles and all internal 
edges. Then by the definition of <ph and V^,, the scaling argument, and the Schwarz 
inequality, 



\{n h Vv-V w R h v,V w (t> h )\ 



^ (qft, Vw(<f>b\e))K e 



E Ob, <l h • n) e 

£ E h WML^(e)\hh\\L^(e) 

^ l|qfelU~(n) E ft (H0o||L°°(jc e 
^ l|qfelU»(Q) E ll^lloA.Jfe 



1100 



»b||L~(e)J 



< 



||qfe|U~(n) 



E II^Ho,*,*. 




</i 2 ||qfelU»(n)ll0A||o,/t. 
Now, by inequalities (|3.2p and (J373TJ) , we have 

||qhlU«»(a) < l|v« - n h Vu|| L oc (0) + || v« - v u ,i?/ l w|| L =c (0) 

<h n+1 \\v\\ wn+ ^ m + h n+1 \lnh\\\Av\\ w ^ m , 

for < n < j. This completes the proof of part (i). 

The proof for part (ii) is similar. One simply needs to replace inequality (|3. 371) 
by (|3.38p in the estimation of ||q?i||L o 



(n)- 



4. Error analysis. The main purpose of this section is to analyze the approx- 
imation error of the weak Galerkin formulation (|2.4[) . For simplicity, in this section, 
we assume that the solution of (|2.4I) satisfies u G £f 3+7 (f2) and w G ff 1+7 (fi), where 
7 > This is not an unreasonable assumption, as we know from p.2|) . the solution 
it can have up to H 4 regularity as long as ft satisfies certain conditions. However, our 
assumption does not include all the possible cases for the biharmonic equation. 

Testing w = —Au with <f>h = {0o, </>&} G Vh and then by using (I3.7[) we have 

(4.1) (w, 4> h )) = (to, O ) = -(V • Vu, 0o) = (n ft Vu, V w cj) h ). 
Similarly, testing —Aw = / with tph ~ {^OiV'f)} G ^o,/i gives 

(4.2) (n h V«;,V to ^) = (/,^o). 
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Comparing (I4.1[) - (|4.2p with the weak Galerkin form (12.41) . one immediately sees that 
there is a consistency error between them. Indeed, since Vh and V ,h are not sub-spaces 
of H 1 ^) and Hq(£1), respectively, the weak Galerkin method is non-conforming. 
Therefore, we would like to first rewrite (|4.1I) - (|4.2I) into a form that is more compatible 
with (j2l| . By using (j3~26| and ((337|) . equations (|LTj) - (l4~2|) can be rewritten as 



(4.3) 
where 



lN h w, <p h )) - (V w R h u, V h (f> h ) = E(w, u, (f) h ), 
(V w N h w,V w ip h ) = {f,i>o), 



E(w, u, 4>h) = lN h w - w, 4>h} + (n^Vu - V w Rhu, V w (j>h)- 



Dchnc e u = Ruu — Uh G Vb.ft, and e w = NhW — Wh G Vh- By subtracting (I4.3P 
from (12.41) . we have 



(4.4) 



{e w ,(f>h} - (V w e u ,V h (f) h ) = E(w,u,4>h) for all (f> h G Vh, 

(V w e w ,V w ^h) = for all ip h 6 Vo,h- 



Notice here (V w e w , V w iph) = docs not necessarily imply e w — 0, since the equation 
only holds for all iph G Vo,h while e w is in Vh- 

Lemma 4.1. The consistency error E(w,u, <p>h) is small in the sense that 
\E(w,u,(f> h )\ < h m \\w\\ m+1 \\(f)h\\o,h + h n+ i | In h\ |Hvr»+».°°(n)||0fc||o,h, 

where ^ < m < j + 1 and < n < j. Moreover, for j > 1, we have the improved 
estimate 

\E(w,u,<j)h)\ < /i" l ||u;|| m+ i||^||o,/ l + /i™ + 5||u|| H /"+2,=o (0) ||^||o,h, 

where h < m < j + 1 and 1 < n < j. 

Proof. The proof is straight forward by using the Schwarz inequality, Lemma 
Remark |3.4[ and Lemma [3~51 □ 



To derive an error estimate from (|4.4[) , let us recall the standard theory for mixed 
finite clement methods. Given two bounded bilinear forms a(-, •) defined on X x X 
and &(•, •) defined on X x M, where X and M are finite dimensional spaces. Denote 
XoClby 

X a = {</> G X : 6(0, ip) = for all t/> G M}. 
Then for all x e ^ and £ G M, 

sup | U |i ,11,11 ^ X X + |m, 

4>eX,^eM \m\x + \\ip\\M 

if and only if 

su p irnr- ~ llxlk, for all x e x , 

rN 0ex o 1101k 

(4-5) , . 

SU P tSt 2 £ IKIk> for a11 S e M - 

II0IU 
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In our formulation, we set X — Vh with norm || • ||q h an d M = Vo t h with norm 
I • |||. Define 

a{x,<t>) = lx,<P)), K&o = -(y w <j>, v„0- 

It is not hard to check that both of these bilinear forms are bounded under the given 
norms. In particular, the boundedness of &(•,•) has been given in ()3.21j) . It is also 
clear that the first inequality in (|4.5[) follows from the definition of a(-, •) and || • \\o,hi 
and the second inequality follows directly from (|3.22[) . Combine the above, we have 
for all x £ Vh, and £ E V ,h, 

(4-6) sup — jjr-iji > \\x\\o,h + III C III - 

<t>ev h ,4evo, h \m\o,h + WWW 

Theorem 4.2. The weak Galerkin formulation \2.1$ for the biharmonic problem 
has the following error estimate: 



\\£wh,h + lll^lll < h m \\w\\ m+1 + h n+ ?\lnh\\\u\\ W n+2,o 



(Q), 



where i < m < j + 1 and < n < j. Moreover, for j > 1, we have the improved 
estimate 

\\£w\\o,h + IMI < h m \\w\\ m+1 + h n+ ^\\u\\ Wn+ 2,^ m , 

where ^ < m < j + 1 and 1 < n < j '. 
Proo/. By (J33J) and (|4~6)) , 

^ ((<Hu, <K)) - (Vw4>h, y w £u) ~ (V W E W , V w lj)h) 

\\£wh,h + IMI < SUp ij-T-ji 

^eVh.^eVb.h II0mo,/i + lllwll 

E(w,u,(j> h ) 

= SUp — — j: . 

<l>h£V h ,i> h ev , h \\<Ph\\o,h + W\VhW\ 
Combining this with Lemma 14. 1[ this completes the proof of the theorem. □ 

Remark 4.1. Assume that the exact solution w and u are sufficiently smooth. It 
follows from the above theorem that the following convergence holds true 



\\ £ w\\o,h + ||kt, 



< 



r O(h*\\nh\) forj = 0, 
0(h j+ ^) for j > 1. 



At this stage, it is standard to use the duality argument and derive an error 
estimation for the L? norm of e u . However, estimating ||e„||o,/i is not an easy task, 
as is similar to the case of Poisson equations. For simplicity, we only consider ||e Ul o||, 
where e u is conveniently expressed as e u — {e u ,o, Define 



£ + An = 0, 



(4.7) 

where -q — and ^ = on dfl. We assume that all internal angles of ft are less than 
126.283696 • • -°. Then, according to ([2~2"]h the solution to P~7]) has H 4 regularity: 

11^12 + Nk<lkf,o||. 
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Furthermore, since such a domain f2 is convex, the Poisson equation with either the 
homogeneous Dirichlet boundary condition or the homogeneous Neumann boundary 
condition has H 2 regularity. 

Clearly, Equation (|4.7[) can be written into the following form: 



(4.8) 



iNht, M - (VwRhV, V w fa) = n, fa) for all fa = {fa, (f> b } G V h , 
(V w Nh£, V w iph) = (e u ,o, i>o) for all iph = {ipo, ^fc} S Vo.fc. 



For simplicity of the notation, denote 

A(N h £,R h r) h ; fa,ip h ) = (N h £, fa} - {V w R h V, V w fa) - (V„,JV fc f, V w ip h ). 

Note that A is a symmetric bilinear form. By setting fa — e w and iph = e u in 
and then subtract these two equations, one get 

||£u,o|| 2 = E(£, T), e w ) - A(N h £, R h r); s w ,e u ) 
(4.9) = E(£ Nh^RhV) 

= E(£,ri,e w )-E(w,u,N h £). 



Here we have used the symmetry of A(-, ■) and Equation (14.4[) . 

The two terms, E(£,r),s w ) and E(w,u, Nh£), in the right-hand side of Equation 
(|4.9[) will be estimated one by one. We start from E(£,n,e w ). By using Lemma ETT1 
it follows that 

(i) When j = 0, 

^ 41Q ^ E(£,r),e w )< (/i||^|| 2 + /i5|ln/i|||7 y || B ,2,=„ (n) J He^Ho,* 
<^ 1/2 |ln/i|(||el|2 + hl|4)||^||o,h. 

(ii) When j > 1, let 5 > be an infinitely small number which ensures the 
Sobolev embedding from W^ 2 {VL) to W 3 -^°°(n). Then 

^ n , E(Z,t),e w ) < ^/i||^|| 2 + /!,3- 5 |ln/i|||7/|| W 3- 5 ,c» ( n)J \\e w \\o,h 

<h(Uh+HU)\\e w \\o, h - 
Next, we give an estimate for E(w, u, Nht;). 

Lemma 4.3. Assume all internal angles oftt are less than 126.283696 • • °, which 
means the biharmonic problem with clamped boundary condition in ft has ft 4 regular- 
ity. Then 

(i) For j = , 

E(w,u,N h 0< (^|k|| m+ i+^||(/-Qo)/||+^ +1 h||„ + i) ||£||a, 

where i < m < 1 and 1/2 < n < 1. 
(ii) For j > 1, 

E(w )U) N h < (h m+1 \\w\\ m+1 + h 2 \\(I - Q )f\\ + h n+1 \\u\\ n+1 ) ||£|| 2 , 
where ^ < m < j + 1 and 1/2 < n < j ' + 1. 
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Proof. By definition, 
(4.12) E(w,u,N h £) = (N h w - w,N h O + (n h \7u - V w R h u,V w N h £). 

First, by the definition of ((•,•)), the Schwarz inequality, Remark 13.31 and [3. 41 we have 
(N h w-w,N h 0) 
= (N w-Q W,N + J2 h{N w - N b w,N £ - N b £) dK 

K£T h 

(4 - 13) < \\N Q w - QoHIHJVbeil + \\N h w - w\\ , h \\N h t - £\\o,h 

'(h m +i\\w\\ m+1 + h 2 \\(I-Q )Aw\\)U\\2 forj=0, ±<m<l 
(fc m+1 |Mlm+i + h 2 \\(I - Q )A«;||)||f||2 for j > 1, i < m < j + 1 ' 



< 



Next, by using inequalities IpTgjl . ff3T27|) . (j3~7| . (j3~4|) . ([3"3Tj) and (pH2|) one after one, 
we get 

(II/, Vu - V w R h u, V w N h £) 
= ((II/, - Ph)Vu, V^JVhf) + (V™(Q/,u - R h u),V w N h £) 
= ((II/, - P ft )Vu, V^iV ft + (V w (Q fe u - flhu),n h V0 

= ((II,, - P A )Vti, V^iV^ - Q/,0) + (( n '> - P >0 Vu > p hV0 - (Qou - Rou, A£) 

<h n+1 \\u\\ n+1 u\\2 + ((n h - j)v«,p & vo + h 2 \\(i - Qo)At»||||e|| a , 

for i < n < j + 1 . The estimation for ((11/, — I)V«, P/, V£) follows the same technique 
used in Inequality (|3.34[) . By the definition of II;, and since |^ = on c?f2, we know 
that (Tib — I)\7u ■ n also vanishes on dfl. Therefore, using the same argument as in 
(|3.34p . one has 

((U h - I)S7u,P h VO < h n+1 \\u\\ n+1 U\\2 + h 2 \\(I - Q )A«||||e|| 2 
for i < n < j + 1. Combining the above gives 

(4.14) (Tl h Vu - V w R h u, V w N h £) < (h n+1 \\u\\ n+1 + h 2 1| (I - Q )Au\\) ||£|| 2 . 

for | < n < j + 1 . 
Notice that 

h 2 \\(I - Q )Au\\ = h 2 \\(I - Q Q )w\\ < h m+2 \\w\\ m for0<m<j + l, 

fc 2 ||(J-Q )AHI=/i 2 ||(/-Oo)/||. 

The lemma follows immediately from (14. (|4. 15[) . □ 

Finally, combining Theorem 14. 2 [ inequalities (|4.9j) . f|4.10f) - (j4.11[) . and Lemma l4~3l 
we get the following L 2 error estimation: 

Theorem 4.4. Assume all internal angles of fl are less than 126.283696 •• ° , 
which means the biharmonic problem with clamped boundary condition in O has H A 
regularity. Then 
(i) For j = 0, 

IKoll < h m+ i\lnh\\\w\\ m+1 + h\hih\ 2 \\u\\ w z,oo (n) 
+ h 2 \\(I-Q )f\\ + h n+1 \\u\\ n+l , 

where \ < m < 1 and i < n < 1. 
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(ii) For j > 1, 

Koll < h m+1 \\w\\ m+1 + h l+ i\\u\\ w i + 2,oo m + h 2 \\(I - Q )f\\ + h n+1 \\u\\ n+1 

where h < m < j + 1, i < n < j + 1 and 1 < I < j '. 
Remark 4.2. //u 7 it) and / are sufficiently smooth, then we get 



\Su,o\\ 



< 



0{h\\nh\ 2 ) forj = 0, 
0(hP + §) forj>l. 



5. Numerical results. In this section, we would like to report some numerical 
results for the weak Galerkin finite element method proposed and analyzed in previ- 
ous sections. Before doing that, let us briefly review some existing results for H 1 -!! 1 
conforming, equal-order finite element discretization of the Ciarlet-Raviart mixed for- 
mulation. As discussed in [5J |3FJ, theoretical error estimates for such schemes are 
indeed sub-optimal due to an effect of inf X(i \\u — XfclU, where Xh is taken from the 
employed H 1 conforming finite clement space. For example, when H 1 -!! 1 conform- 
ing quadratic elements are used to approximate both u and w, the error satisfies 
||u - Uhh + \\w -Wh\\ < inf Xh \\u - Xhh + inf xh \\ w ~ Xh\\ < 0(h), while intuitively, 
one may expect \\w — Wh\\ to have an 0(h 2 ) convergence. By using the L°° argument, 
Scholz [48) was able to improve the convergence rate of L 2 norm for w by hh , and it 
is known that this theoretical result is indeed sharp. For the weak Galerkin approx- 
imation, from the discussing in the previous sections, clearly we are facing the same 
issue. 

However, numerous numerical experiments have illustrated that H 1 -!! 1 conform- 
ing, equal-order Ciarlet-Raviart mixed finite element approximation often demon- 
strates convergence rates better than the theoretical prediction. Indeed, this has been 
partly explained theoretically in [49], in which the author proved that optimal order 
of convergence rates can be recovered in certain fixed subdomains of fi, when equal 
order H 1 conforming elements are used. We point out that similar phenomena have 
been observed in the numerical experiments using weak Galerkin discretization. This 
means that numerical results are often better than theoretical predictions. 

Another issue in the implementation of the weak Galerkin finite element method 
is the treatment of non-homogeneous boundary data 

u = gi on dtt, 

du an 
— = g 2 on di.1. 

Clearly, both boundary conditions are imposed on u, and u = gi is the essential 
boundary condition while ^ = gi is the natural boundary condition. To impose the 
natural boundary condition, we shall modify the first equation of (|2.4[) into 

((Wh, 4>h} - (VwUh, V w (t>h) = -(52, 4>b)dQ.- 

The essential boundary condition should be enforced by taking the L 2 projection of 
the corresponding boundary data. 

Consider three test problems defined on 57 — [0, 1] x [0, 1] with exact solutions 

ui =x 2 {l-x) 2 y 2 {l~y) 2 , 

U2 = sin(27rx) sin(27ri/) and u% = sin(27rx + — ) sin(27ry + — ), 
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respectively. The reason for choosing these three exact solutions is that they have the 
following type of boundary conditions 

ui\dn = 

u 2 \dn = 

U3\dn ^ 

This allows us to test the effect of different boundary data on convergence rates. 
Although the theoretical error estimates are given for e u = RhU — Uh and e w = 
NhW — Wh, it is clear that they have at least the same order as e u — QhU — 
and e w — Nh.w — Wh, provided that the exact solution is smooth enough. Thus for 
convenience, we only compute different norms for e u and e w , instead of for e u and s w . 

The tests are performed using an unstructured triangular initial mesh, with char- 
acteristic mesh size 0.1. The initial mesh is then refined by dividing every triangle 
into four sub-triangles, to generate a sequence of nested meshes with various mesh size 
h. All discretization schemes are formulated by using the lowest order weak Galerkin 
element, with j = 0. For simplicity of notation, for any v G Vh, denote 

Nil = ( £ h\\v b f a 
\Ker h 

The results for test problems with exact solutions u±, u 2 and «3, are reported 
in Table 15.11 15.21 and 15.31 respectively. The results indicate that u always achieves 
an optimal order of convergence, while the convergence for w varies with different 
boundary conditions. It should be pointed out that both of them have outperformed 
the convergence as predicted by theory. 



Table 5.1 

Numerical results for the test problem with exact solution u\ and lowest order of WG elements. 



h 


||V lu e u || 


IK.oll 


\\ e u,b\\ 


| ^ ' W&W || 


l|e»,o|| 


\\ e w,b\\ 


0.1 


1.33e-03 


2.40e-04 


4.59e-04 


5.66e-02 


2.96e-03 


6.91e-03 


0.05 


4.69e-04 


6.18c-05 


1.17e-04 


2.80e-02 


9.14c-04 


1.99e-03 


0.025 


2.00e-04 


1.55e-05 


2.97e-05 


1.60e-02 


2.64e-04 


5.70e-04 


0.0125 


9.56e-05 


3.90e-06 


7.44e-06 


1.21e-02 


8.33e-05 


1.89e-04 


0.00625 


4.72e-05 


9.77e-07 


1.86e-06 


1.13e-02 


3.26e-05 


7.91e-05 


Asym. Order 
0(h k ), k = 


1.1930 


1.9876 


1.9877 


0.5864 


1.6461 


1.6298 



dn 

du 2 
dn 
du 3 
dn 



= 0, 



Ml 



mi 



= 0. 




Our final example is a case where the exact solution has a low regularity in the 
domain ft = [0, 1] x [0, 1]. More precisely, the exact solution is given by 

Q/2 I . 30 . 0\ 

114 = r ' sin 3 sin - , 

V 2 2) 

where (r, 9) are the polar coordinates. It is easy to check that u S H 2 5 . The errors 
for weak Galerkin finite element approximations are reported in Table 15.41 Here, u 
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Table 5.2 

Numerical results for the test problem with exact solution U2 and lowest order of WG elements. 



h 


||V w e u 


l|e«,o|| 


l|e«,6|| 


\\V w e w \\ 


l|e«,,o|| 


\\e w ,b\\ 


0.1 


9.58c-01 


8.66c-02 


1.65c-01 


4.39c+01 


6.09c-01 


2.01c+00 


0.05 


3.34c-01 


2.18o-02 


4.14c-02 


2.32c+01 


2.78c-01 


7.19c-01 


0.025 


1.43c-01 


5.47e-03 


1.03c-02 


1.37c+01 


1.15c-01 


2.81c-01 


0.0125 


6.81c-02 


1.37e-03 


2.59c-03 


1.02c+01 


5.12c-02 


1.26e-01 


0.00625 


3.36c-02 


3.42e-04 


6.49c-04 


9.33c+00 


2.45c-02 


6.12c-02 


Asym. Order 
0{h k ), k = 


1.1958 


1.9958 


1.9975 


0.5649 


1.1709 


1.2587 



Table 5.3 

Numerical results for the test problem with exact solution U3 and lowest order of WG elements. 



h 


\\V w e u \\ 


IK,o|| 


IK&II 


||V„,e w || 


IK,o|| 


\\e w ,b\\ 


0.1 


8.23c-01 


1.18c-01 


2.27c-01 


5.61c+01 


4.25c+00 


9.42c+00 


0.05 


3.07c-01 


3.18c-02 


6.09e-02 


2.43c+01 


1.24c+00 


2.58c+00 


0.025 


1.35c-01 


8.13c-03 


1.55c-02 


1.13c+01 


3.28e-01 


6.61e-01 


0.0125 


6.49e-02 


2.04c-03 


3.90e-03 


5.58c+00 


8.42e-02 


1.67e-01 


0.00625 


3.21c-02 


5.11c-04 


9.78e-04 


2.77c+00 


2.14c-02 


4.21c-02 


Asym. Order 
0(h k ), k = 


1.1599 


1.9679 


1.9682 


1.0801 


1.9157 


1.9558 



still achieves an optimal order of convergence, while the convergence rates for w is 
restricted by the fact that w e H 5 . All the results are in consistency with the theory 
established in this article. 



Table 5.4 

Numerical results for the test problem with exact solution «4 and lowest order of WG elements. 



h 


\\V w e u \\ 


l|e«,o|| 


l|e«,b|| 


llV^e^H 


IK,o|| 


\\e w ,b\\ 


0.1 


3.73e-02 


9.44e-04 


2.15c-03 


2.88e+01 


4.05c-01 


1.78c+00 


0.05 


1.87c-02 


2.55e-04 


5.73c-04 


4.08c+01 


2.86c-01 


1.26c+00 


0.025 


9.37c-03 


6.60c-05 


1.46e-04 


5.77c+01 


2.02e-01 


8.91c-01 


0.0125 


4.68c-03 


1.67c-05 


3.69e-05 


8.16c+01 


1.42e-01 


6.30e-01 


0.00625 


2.34c-03 


4.19e-06 


9.24c-06 


1.15c+02 


1.01e-01 


4.45c-01 


Asym. Order 
0(h k ), k = 


0.9984 


1.9567 


1.9690 


-0.4998 


0.5008 


0.5000 
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